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An efficient density matrix renormalization group (DMRG) algorithm is presented for the Bethe 
lattice with connectivity Z = 3 and antiferromagnetic exchange between nearest neighbor spins 
s = 1/2 or 1 sites in successive generations g. The algorithm is accurate for s = 1 sites. The ground 
states are magnetic with spin S{g) — 2®s, staggered magnetization that persists for large g > 20 and 
short-range spin correlation functions that decrease exponentially. A finite energy gap to 5* > S{g) 
leads to a magnetization plateau in the extended lattice. Closely similar DMRG results for s = 1/2 
and 1 are interpreted in terms of an analytical three-site model. 
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I. INTRODUCTION 

Dendrimers have been extensively explored over the 
past two decades, both experimentally and theoretically. 
Dendrimer are macromolecules with repetitive branches 
attached to a central core or focal point ^ ■ Three pa- 
rameters (b,Z,g) characterize the dendrimer lattice model 
in FigHJ the number of branches b attached to the fo- 
cal point, the connectivity Z of a site and the number 
of generations g. Many natural systems are dendrimers 
and molecular engineering can produce custom den- 
drimers with constituents that range from carbon-based 
molecules [1] to organo- metallic compounds [3, Po- 
tential applications include drug delivery for chemother- 
apy 3 , gene therapy , magnetic resonance imaging as 
contrast agents based on superparamagnetism [7], and 
molecular recognition Q. 



Bethe lattices (BL) or Cayley trees are theoretical 
models of dendrimers and attractive approximations of 
solid-state systems. We consider in this paper super- 
paramagnetic dendrimers based on sites with unpaired 
spins. BL's are fundamentally one-dimensional (ID) sys- 
tems, without any closed loops, that can be divided by 
cutting any bond. Their distinctive feature is exponential 
growth with g. Half of the sites for the BL with Z = 3 
in Fig. [T] are on the surface in generation g. BL's are 
bipartite: all nearest neighbor (NN) bonds in FigH] are 
between sites in even and odd g that define two sublat- 
tices with different numbers of sites. The BL has open 
boundary conditions, a focal point at g — and an arbi- 
trarily large boundary g. General theorems apply for the 
spectrum free electrons ^ with NN transfer t in a BL or 
for spins [lOj with NN Heisenberg exchange J > 0. Three 
parameters generate a rich variety of BL(b, Z, g). In dy- 
namical mean field theory , the BL density of states 




FIG. 1: One arm of a Bethe lattice with b = Z = 3 and 
g = 5. The other two arms are shown to generation g = 1; g 
also represents distance from the focal point. 



with infinite connectivity is used as an initial guess for 
the density of states in higher dimensions. BLs are mod- 
els for strongly correlated systems [l^, alloys and 
disordered systems ■ The present study addresses an- 
tiferromagnetic (AF) Heisenberg exchange J > in the 
BL with b = Z = 3 
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The sum is over all NN for either s = 1/2 or 

s — 1 sites. BL-Ising models limited to s| interactions 
have been used to study finite spin glasses [l^. The 
electronic properties of correlated BL models still pose 
many challenges. 

Recent advances in numerical techniques and com- 
putational resources have been applied to BL models. 
Methods include exact diagonalization (ED), quantum 
Monte Carlo (QCM) [H] and Density Matrix Renor- 
malization Group (DMRG) DMRG is particularly 

well suited for ID systems such as Hubbard or extended 
Hubbard models, t — J models and Heisenberg or related 
spin models. DMRG yields accurate properties for the 
ground state (gs) or low-energy excited states 18-21]. 
The DMRG challenge for a BL is the large number of 
surface sites in Fig. [TJ Otsuka [l^] applied DMRG to 
the BL, Eq. [TJ with s = 1/2 sites and axially anisotropic 
NN exchange (XXZ model). Friedman [23| presented 
another DMRG algorithm for Eq. [T] with s = 1/2. 
Lepetit et al. [l^l used essentially the same algorithm to 
treat the Hubbard model version of Eq. [T] that reduces 
to J = At^ /U for s = 1/2 when the on-site repulsion 
U is large compared NN electron transfer t. They also 
solved analytically the Hiickel or tight-binding model 
with U = 0. 




FIG. 2; Schematic representation of BL growth from g = 1 
at step 1 to g = p at step p. The blocks A = B = C at step g 
are given by Eq. (2] The focal point D is the new site added 
at each step. 



DMRG is a truncation procedure in which insignifi- 
cant degrees of freedom of the system block, the right 
or left block in a chain, are discarded at each step with 
increasing system size. In ID chains, the superblock 
consists of two blocks with dimension ni and two sites 
with p degrees of freedom. The superblock dimension 
is ni? X p^, with p = 2s + 1 in spin systems. The 
superblock of BL(3, 3, g) has three blocks and hence goes 
as rrv^ . More branches h increases the computational 
requirements. Otsuka [22| used four blocks and two new 
sites for BL(3, 3, g) with p = 2 {s = 1/2). His su per block 
increased as x p^ . The Friedman algorithm [23| with 
h — Z and four new sites yields a superblock dimension 
of X p^. The DMRG algorithm in Section II has a 
superblock dimension m? x p or, more generally, x p, 
that makes s = 1 sites accessible. 

The paper is organized as follows. Section II presents 
and tests the new algorithm. Section HI reports results 
for the BL in Fig. [Hand Eq. [T]up to g = 11 for s = 1/2 
and 1 sites, and up to g > 20 for (sg) at the focal point. 
We obtain the gs energy per site, the energy gap A that 
governs the gs magnetization and gs expectation values of 

for n < g. We find spin correlation functions between 
the focal point and other sites and their convergence with 
increasing g. Section IV relates the DMRG results to a 
simple analytical model with localized states, to previous 
DMRG studies and to the question of long-range order 
in the infinite BL. 



II. DMRG ALGORITHM 

In this section, we present a new DMRG algorithm for 
the BL in Fig. [T]with b — Z ^ 3. The principal change 
is how the lattice is grown. The total number of sites in 
BL(3,3,g) is 

Nrig) = 1 + 3(25 _ i) (2) 

The first step in Fig. [51 contains four sites. Sites A, B 
and C are blocks each of whose size is 2^"^ at generation 
g. The focal point D is the new site added at each step. 
As shown in Fig. [2l each block A = B = C contains three 
sites in the second step for g — 2. Growth at step g = 2, 
3, ... is schematically represented as 




Dg is connected to the focal point D at step g + 1- The 
next step gives blocks of 7 sites at g = 3, and so on up to 
the BL with NT{g) sites. The procedure for BL(3,3,g) 
holds for other b and Z. 

The infinite DMRG algorithm for the BL proceeds 
along largely standard lines T?]: 1. Start with the su- 
perblock matrix of four sites and find the eigenvalues of 
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TABLE I: The difference 8X between ED and DMRG with 
increasing m for BL(3,3,3) with s=l/2 sites, where X is the 
gs energy per site eo, excitation energy A, the boundary spin 
(s|) and spin correlation (sfsa) 



m 


5e X 10^ 


5A X 10" 


5{sl) X 10" 


5{slsl) X 10» 


10 


2670.8 


42.593 


89.48 


-1.193 


20 


560.08 


41.441 


0.030 


-0.722 


30 


552.99 


41.212 


0.031 


-0.801 


50 


0.001 


0.0008 


0.003 


-0.0001 



H. 2. Use the eigenvectors of the superblock to construct 
the density matrix of the new blocks ^g+i, initially for 
(7 = 1. Keep the eigenvectors of the m largest eigenvalues, 
with m chosen as discussed below. The density matrix 
dimension is w? x p, where m and p refer to the block 
and degrees of freedom of the new site. Full diagonal- 
ization of the density matrix is carried out separately for 
large m in sectors with different total . 3. Renormal- 
ize the Hamiltonian of the new blocks and the operators 
that are necessary for the next step. These steps follow 
conventional DMRG [l^. 4. Construct the next {g + 
1) superblock from the three renormalized blocks Ag^i 
and the new site. Diagonalize the matrix and retain the 
m lowest eigenvalues and eigenvectors. Repeat steps 2-4 
until the desired system size is reached. 

The superblock dimension of x p makes possi- 
ble larger m that increases the accuracy and larger 
p = 2s -|- 1. We can use to = 60 without much computa- 
tional effort and find 10^^'^ or less for the weight of the 
discarded eigenvalues of blocks. DMRG is a variational 
method. Energies and correlation functions for given size 
g converge better for finite DMRG [l7[ . We followed the 
standard approach of sweeping back and forth through 
different blocks. Care has to be taken in designing the fi- 
nite DMRG algorithm due to the BL's complex structure. 

We constructed the density matrix with equal weight 
for the lowest two eigenstates. As a first test of accuracy, 
we performed ED on BL(3,3,3) with 22 sites s=l/2 and 
on BL(3,3,2) with 10 sites s=l. There are 2^^ ^nd 3^° 
spin states, respectively. DMRG results with increasing 
TO must eventually converge to ED. The evolution of the 
gs energy per site, 5tQ = to{m) — eQ{ED), with m is 
shown in Table H] for the 22 sites systems. Also shown 
are the evolution of d{s^) and the spin correlation (5(sgs|) 
between the focal point and the boundary. DMRG with 
TO = 50 is quantitative here. The Friedman algorithm 
whose super block increases as x p'^ is limited to to « 
30 and returns ^ Seo = 2.2 x 10"^ with to = 29 for 
infinite DMRG, which is comparable to m = 10 in Table 
m DMRG with increasing to also agrees quantitatively 
with ED for s=l sites. 
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FIG. 3: Ground state energy Eo of spins 1/2 and 1 in a large 
BL(3,3,10) as a function of m relative to m = 65. The inset 
shows the m-dependence of the energy gap A to the ground 
state with = ^(lO) + 1 



We chose to = 60 on the basis of BL calculations 
for large systems with g — 10. Figure [3] shows the 
evolution of eo(TO)/eo(65) with m for s = 1/2 and 1 sites. 
Although s=l converges an order of magnitude more 
slowly, TO > 40 is adequate in either case. Preliminary 
results indicate still slower convergence for s=3/2 sites. 
The inset of Fig. |3] shows the slower, and not the 
monotonic, evolution of energy gap A(to)/A(65). This 
is not unexpected since A is the difference between two 
extensive quantities. 

We followed the first excited state in = S{g) with 
similar results. We kept to = 60 and did 8 sweeps of finite 
DMRG for the results in Section III. We estimate that the 
gs energy, spin densities and spin correlation functions 
are accurate to 4-5 decimal places in larger systems, while 
energy gaps are accurate to 2-3 places. 



III. RESULTS FOR BL WITH S = 1/2 AND 1 

In general, a BL of Nt sites has Nt — 1 bonds, since 
only the focal point is not connected to a site with lower 
g. The gs energy per bond of BL(3,3,.g) is 



eo(g) - Eo{g)/{NT{g) - 1) 



(4) 



where Eq is the gs energy and Nt is given in Eq. [2] 
Fig. [5] shows eo{g) for s = 1/2 and 1 sites up to g = 
10. There is very little size dependence. The extended 
s = 1/2 and 1 systems have eo=-0.39384 and -1.2795, 
respectively. 

The gs has total z component of spin S^'{g) = 2^s. 
Eq. [T] conserves S and either a half-filled band or 
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FIG. 4: Ground state energy eo per bond of the BL, Eq. [T] 
with s = 1/2 or 1 sites, g < 11 generations and Nrig) sites. 
Linear extrapolation gives eo of the extended BL. 



Heisenberg exchange yields S{g) = 2^3. The gap A{g) 
to the gs in the = S{g) + 1 sector governs the gs 
magnetization. The evolution of A((7) with Nxig) is 
shown in Fig. [SJ The substantial gap of the extended 
system is discussed in Section IV. We also found doubly 
degenerate excitations A'{g), comparable to A(g), in the 
S^{g) sector for both s = 1/2 and 1. The C3 symmetry 
of the BL leads naturally to E states. Since we compute 
rather than S, the gs automatically appears also in 
sectors with < S^{g). The second and third excited 
states for = S^{g) — 1 decrease with Nrig) and vanish 
in the extended system within our numerical accuracy. 

The gs expectation values of (sp in generation p 
are listed in Table M for the BL with g = 10. The g 
— 9 and 11 results are almost the same. The gs has 
long-range order (LRO) that corresponds to staggered 
magnetization in successive generations. The largest 
(Sp) is at the boundary, p = g. The smallest magnitude 
is a,t p — g — 1 next to the boundary, and |(s^)| near the 
focal point become equal for large g. The convergence of 
|(sg)| to 0.348 is shown in Fig. IHlup to g = 26, a huge BL 
of 3 X 2^^ sites, and agrees with the previous estimate of 
0.35 [24]. The s = 1 limit up the g = 24 is |(sg)| = 0.83, 
and the staggered magnetization of both gs becomes 
constant near the focal point for g > 10. By contrast, 
a half- filled BL of free electron has (sg) = 2/{g + 1) for 
odd g and (sq) = for even g or in the extended system. 
Table mi also lists the gs expectation values of (sp_i ■ Sp) 
in successive generations of BL(3,3,10) with s = 1/2 
and 1 sites. The g = 9 and 11 values are similar. 
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FIG. 5: Size dependence of the energy gap A to 5^ = S(g) 
-I- 1 for the BL, Eq. [T] with s = 1/2 or 1 sites, g < 11 and 
Nrig) sites. 



Variations of (sp-i • Sp) at g and g — 1 are reduced near 
the focal point. The BL-Ising model has (Sp) = (— 

1 



and 



Table |TT] shows that s 



IS 



closer to the Ising model than s — 1/2, where quantum 
fluctuations are larger. 



We define radial spin correlation functions as 
C{r) = {-in{sh-s;)-{s^,){s^^)) 



(5) 



with r = 1,2. ..5. The mean-field contribution is ex- 
plicitly excluded. As seen in Fig. [3 C(r) decreases as 
exp(— ar) with a = 0.80 for both s = 1/2 and 1. Ac- 
curate DMRG makes it possible to compute small C(r) 
up to 7' « 10. The inset shows the related correlation 
function C^(r) with (sqS^) instead of the dot product, 

TABLE II: Ground state expectation values of (Sp) and (sp-i- 
Sp) for s=l/2 and 1 site in BL(3,3,10). The focal point and 
boundary are p = and p = 10 respectively. 



p 




s=l/2 




s=l 




(4) 


(•^p— 1 ' Sp^ 









0.347 




0.829 




1 


-0.345 


-0.359 


-0.828 


-1.214 


2 


0.347 


-0.361 


0.828 


-1.214 


3 


-0.342 


-0.357 


-0.826 


-1.213 


4 


0.348 


-0.364 


0.826 


-1.215 


5 


-0.332 


-0.350 


-0.817 


-1.208 


6 


0.348 


-0.371 


0.820 


-1.218 


7 


-0.307 


-0.334 


-0.789 


-1.191 


8 


0.356 


-0.392 


0.816 


-1.238 


9 


-0.249 


-0.289 


-0.691 


-1.104 


10 


0.393 


-0.456 


0.872 


-1.397 
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FIG. 6: Magnitude of {sg) at the focal point of a BL up to g 
= 26 generations for s = 1/2 and up to y = 24 for s = 1 sites. 
Even and odd g form separate series that merge at large g. 



which decreases even faster. C(r) is almost completely 
due to transverse spin components. Since the number 
of boundary sites goes as n{g) — 1.5exp(gln2), we have 
n{g)C{g) -> for large g. 

The gs degeneracy is lifted when an applied magnetic 
field h is added to Eq. [T] The lowest Zeeman level goes 
as —hS{g). With increasing h, the gs becomes the lowest 
Zeeman level of a state with > S{g), an excited state 
at /i = 0. When S*^ = S{g) + 1, the crossover field is 
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FIG. 7: Spin correlation functions C(r) in Eq. [5] between 
the focal point and generation r < g in BL's with g = 9, 10, 
11 and s = 1/2 or 1 sites. The inset shows the z component, 
C^r). 
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FIG. 8: Exact magnetization M per site of BL(3,3,3) with s 
= 1/2 sites as a function of applied field h. The g = 4 and 5 
magnetization is DMRG up to finite h. The extended BL has 
an M = 1/6 plateau up to /i = A, the gap shown in Fig. [S] 



related to the zero-field energy in the two sectors. 



h^Aig)^[EoiSig) + l)~Eo]/J. 



(6) 



The first crossover may be to a state with higher S^, 
when Eq. [6] has a multiple of h. The gs magnetization 
per s = 1/2 site of BL(3,3,(?) is shown in Fig. [5]for 5 = 3, 
4 and 5. The first jump is to = S{g) + 1. Complete 
alignment at large h leads to M = 1/2 in reduced units. 
The extended BL has a magnetic gs with M = 1/6 at /i 
~ and an initial increase at A marked with an arrow 
in Fig. m 



IV. DISCUSSION 

DMRG results for H{g) in Eq. [T] are similar for 
BLs with s — 1/2 and 1 sites in Fig. [TJ The gs has 
S{g) = 2^3 as expected on general grounds. The larger 
sublattice has Na{2) = 2P^^ — 1 sites and the smaller 
has NB{g) = 29-1 sites. The difference NA^g) - Nsig) 
is 2^. As seen in Table ITTl the gs has staggered magneti- 
zation with LRO and AF spin correlations in successive 
generations. Spin correlation functions C(r) in Fig. [8] 
decay rapidly and exponentially. Magnetization has a 
substantial gap A.{g) that remains finite in the extended 
s = 1/2 or 1 system. 

We interpret these results in terms of a simple analyti- 
cal approximation. We partition H{g) with J = 1 in Eq. 
m^^s H = Ho + V, 
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Ho = NB{g)h{3) 
his) = si ■ (si + si,) (7) 

The trimer h{3) is elementary and Hq accounts for ex- 
actly 2/3 of the exchanges. For example, the outermost 
two generations of BL(3,3,5) in Fig. [T] contain eight 
trimers per arm. Each s — 1/2 trimer has a gs with 
s = 1/2 at eo — —1, another doublet at e = and a 
quartet, s — 3/2, at e = 1/2. The gs of an s = 1 trimer 
is Co = —3. The perturbation V contains an extra site q 
in the larger sublattice and remaining exchanges. Each 
of N^iq) choices of q uniquely defines the sites of Nsig) 
trimers. The middle and end sites are necessarily in 
the smaller and larger sublattice, respectively, and all 
exchanges in V are between a middle and end site. There 
is only one kind of trimer-trimer or trimer-g interaction. 

The gs energy of Hq for s = 1/2 and 1 sites is —Nsig) 
and —SNsig), respectively. Spin correlation functions of 
adjacent sites for s = 1/2 and 1 are 

(si -5-2) = -1/2, -3/2. (8) 

The spin densities of the central and terminal sites for s 
= 1/2 are 

2(s^) --1/3 , 2(sJ) =2/3, (9) 

The corresponding spin densities for s = 1 are —1 in the 
middle and 3/2 at the ends. The gs of Hq is a product 
over trimers and the spin at q, 

= \a,) n IV',) (10) 

First-order perturbation theory in V lifts the degeneracy 
between q in generations g and < g. 

Site 5 is a 1/Nt correction for large g. The energy per 
site of the infinite BL is, to first order in V, 

(v|/(")|i7|v|/(0))/37V5(5) = eo(s)/3+(.D(^2)/3(ll) 

Parallel spins give lower energy due to the opposite 
signs of spin densities. The s = 1/2 result is -19/54 = 
0.35185 per site in first order and -0.39082 in second 
order in V, very close the DMRG result of 0.39385 in 
Fig. H The s = 1 result is -9/8 = 1.125 per site in 
first order and -1.3152 in second order, slightly below 
-1.2796 in DMRG. The variational theorem holds for the 
energy in first order, but not in second order. The first 
order energy is lowest for parallel spins of Nb (<?) trimers 
and site q that properly gives = S{q)s. Moreover, 
^''''•'(g, (7) immediately rationalizes low-energy spin flips 



in states with S < S{g) for either s = 1/2 or 1 sites. 

Site (7 is a NN of one middle site when q is in 
generation g and of three middle sites otherwise. The 
second term of Eq. [12] is more negative for g < g — 2 
than for q = g hy 2(s|)(s — (sf)) for s = 1/2 or 1. 
To first order in V, site q in not on the boundary and 
g, g — 1 in Fig. [T] are trimers with end sites g. Hq 
accounts for all exchanges between generations g and 
g — 1 while V contains all exchanges between g — 1 and 
g — 2. The trimer approximation leads to (sp = 1/3 
and (Sg_i) = —1/6, compared to 0.40 and -0.25 for 
DMRG for s = 1/2 sites in Table HIl The correlation 
functions are (sg-i • Sg) = —1/2 compared to -0.44 for 
DMRG. Trimers have reduced (sg_2 ■ Sg-i) = —1/18, 
well below the DMRG result of -0.30 but consistent with 
reduced correlation in 5 — 1 and g — 2. The s = 1 BL 
has (Sg) = 3/4 and (Sg_i) = —1/2 in the trimer approx- 
imation and 0.87, -0.69 in DMRG (Table 2). Trimers 
have (sg_i-Sg) = -3/2 for s=l while DMRG gives -1.397. 

A trimer of s = 1/2 sites must be excited to a quartet 
state s = 3/2 with excitation energy 3/2 under Hq to ob- 
tain = S{g) + 1. One of the iVs (5) trimers in ^'('"(q,g) 
is changed to |0) = \aaa). A normalized function with 
a quartet is 

Mg)) = iNsig))-'^' \M\^'^"\q,g))/\^rn) 

7n—l 

(12) 

The quartet is delocalized over the BL by the s+s~ terms 
of V. To first order in V, the excitation energy to 5 = 
S{g) + 1 is 

($|iJ|$) - (*(°)|iJ|«'(")) = A(i)=1.0 (13) 

Delocalization lowers the energy by 2/3 while the 
diagonal s^s^ contribution raises the energy by 1/6. The 
net effect is to lower the excitation from 3/2 to 1 for 
s = 1/2, somewhat above A = 0.74 for the extended BL 
in Fig. [5] 

The sharp distinction between NN exchanges in Hq 
and V for the outermost three generations is lost in 
the interior. Hq contains trimers that spans three 
generations when q < g — 2. When q ^ 0, the central 
site in Fig. [T] is the middle site of a trimer for even 
g and the end site for odd g. Although trimers imply 
intermediate spin density and spin correlations near the 
focal point, there are variations between even and odd g 
in contrast to identical (sg) in Fig. [SI Similarly, C(r) 
in Eq. [5] for 5'(°^((7,f/) is strictly limited to r = 1 or 2 
since trimers span at most span three generations. The 
function "^''■^\q,g) is localized, more localized than the 
DMRG gs, but it rationalizes DMRG results in some 
detail. 
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DMRG results for NN exchange J in Eq. [T] for s = 
1/2 and 1 sites in BL(3,3,g) are closely similar, in sharp 
contrast to the fundamentally different behavior of ID 
chains of s = 1/2 and 1 sites [2^. The s = 1 chain with 
additional NN terms J{si ■ Sj)'^/3 in Eq. [T]is a valence 
bond solid (VBS) with rigorously known gs properties 
[Hi. The VBS on BL(3,3,g) has s = 3/2 site, yet another 
added term to Eq. [TJ and (sg) = and hence no LRO 
in the extended system [2^. Both the spin and Hamil- 
tonian of the VBS are different, and no DMRG has been 
performed on that system. For Eq. 1, DMRG indicates a 
gs with staggered magnetization and finite (sg) in Fig. [5] 
at the focal point for either s = 1/2 or 1. The gs has LRO 
and short-range spin correlations C(r) in Eq. [5] that de- 
crease exponentially and more rapidly than the number 
of boundary sites. 

V. CONCLUSIONS 

The DMRG algorithm in Section II makes it possible 
to treat the BL in Fig. [T]with s = 1 sites. It improves 
the accuracy for s = 1/2 sites and yields low-energy 



excitations. All DMRG results for either s=l/2 or 1 
sites be understood qualitatively in terms of a trimer 
model for BL(3,3,g). 

Preliminary results for BL(3,3,g) with s = 3/2 sites 
are satisfactory. The addition of a focal point with p 
degrees of freedom at each step can be used for BL's with 
more arms 6 > 3 or higher connectivity Z > 3, although 
with steep increases of computational resources. The 
bottleneck is the dimension m^~^ x p of the block whose 
density matrix is constructed at each step; m ~ 70 for 
Z — 3 becomes m « 17 for Z = 4, which has limited 
accuracy. A comparable reduction to m w 20 — 30 
limited previous algorithms [22l - [23 | for BL(3,3,g) to 
s=l/2 sites. The dimension x p of the superblock is 
less serious because only a few eigenvalues are required 
at each step. 
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